Streammetabolizer output for BWL- 30 minute data

Background

This is the lower reach of the larger turbulent mountain stream, Blackwood Creek in CA, USA.

Identifying reasonable modeled estimates of K600

We then ran the normal stream metabolizer model: b_Kb_oipi_tr_plrckm.stan to get modeled K600 to see if we could resolve the negative correlation between ER and K600. Priors on K600_lnQ_nodes_meanlog were set as 5 bins based on rnorm(1000, mean = logQ_mean, sd = logQ_sd) centered around the mean and logQ values 1-2 sd away from the mean.

Load different model segments with normal prior for gas exchange.

We chose segments of time where we believe GPP occurred and was greater than 0. These chunks of time are from a previous model where we binned flow and incorporated measured and estimated K600 priors from gas exchange measurements a the reach.

Metabolism with CIs for the full time series

This the raw model output. It looks okay aside from some small parts of 2023. Where GPP is in blue and ER is in orange, and the black points represent NEP.

Here is the run configuration for full model:

mm_name(type = 'bayes', pool_K600 = "binned", err_obs_iid = TRUE, err_proc_iid = TRUE, ode_method = "trapezoid", deficit_src = 'DO_mod', engine = 'stan')

Fitting priors:

K600_lnQ_nodes_meanlog = log(22) Where 16 was the mean value from observed measurements and normal pooled modeled, K600_lnQ_nodes_sdlog = 1.04 bayes_specs_new$K600_lnQ_nodes_centers <- log_bins was from prior_samples <- rnorm(1000, mean = logQ_mean, sd = logQ_sd)

Checking full model convergence

Make sure the chains converged; all r-hat values were well below 1.05 (the red line) for GPP, ER, and K600. The blue lines are the mean for each parameter.

** Some poor convergence in 2023 for K600

Looks like that weird 2023 time period corresponds to bad rhats for all parameters.

k600.rhat <- mean(na.omit(met.full$K600_daily_Rhat))
k600.rhat

## [1] 1.032055

GPP.rhat <- mean(na.omit(met.full$GPP_Rhat))
GPP.rhat

## [1] 1.00234

ER.rhat <- mean(na.omit(met.full$ER_Rhat))
ER.rhat

## [1] 1.009444

##       date                lab                 rmse               sd          
##  Min.   :2021-04-29   Length:1754        Min.   :0.00944   Min.   :0.009769  
##  1st Qu.:2022-02-04   Class :character   1st Qu.:0.04321   1st Qu.:0.289368  
##  Median :2022-10-12   Mode  :character   Median :0.07265   Median :0.422422  
##  Mean   :2022-12-13                      Mean   :0.10596   Mean   :0.402670  
##  3rd Qu.:2023-10-17                      3rd Qu.:0.15519   3rd Qu.:0.512226  
##  Max.   :2024-08-16                      Max.   :0.69651   Max.   :0.992800  
##                                          NA's   :208                         
##       min              max             range            nrmse        
##  Min.   : 5.339   Min.   : 7.060   Min.   :0.0295   Min.   :0.01029  
##  1st Qu.: 6.971   1st Qu.: 8.407   1st Qu.:0.8977   1st Qu.:0.03958  
##  Median : 8.168   Median : 9.492   Median :1.3170   Median :0.06050  
##  Mean   : 8.067   Mean   : 9.331   Mean   :1.2642   Mean   :0.07767  
##  3rd Qu.: 9.129   3rd Qu.:10.253   3rd Qu.:1.5753   3rd Qu.:0.10953  
##  Max.   :10.578   Max.   :12.078   Max.   :3.8188   Max.   :0.26333  
##                                                     NA's   :208      
##       minT             maxT            rangeT       
##  Min.   : 1.003   Min.   : 1.152   Min.   : 0.0575  
##  1st Qu.: 1.993   1st Qu.: 6.627   1st Qu.: 3.9478  
##  Median : 4.636   Median :10.920   Median : 6.0677  
##  Mean   : 6.165   Mean   :12.052   Mean   : 5.8870  
##  3rd Qu.: 9.931   3rd Qu.:17.508   3rd Qu.: 8.0717  
##  Max.   :16.905   Max.   :26.169   Max.   :10.9748  
## 

Looking at the K600 and ER relationships

Here is the run configuration:

bayes_specs_new

bayes_name_new <- mm_name(type = 'bayes', pool_K600 = "normal", err_obs_iid = TRUE, err_proc_iid = TRUE, ode_method = "trapezoid", deficit_src = 'DO_mod', engine = 'stan')

Where each “lab” segment was run as individual model with the K600_lnQ_nodes_meanlog adjusted to match streamline in cms during that time.

K-ER, K-Q, general K600 trends

There is a strong negative correlation between ER and K600 (-0.96). GPP and K600 are slightly less correlated (0.288).

met.clean <- met.df %>%
  filter(GPP_daily_Rhat<1.1)%>%
  filter(GPP_97.5pct>0)%>%
  filter(ER_daily_Rhat<1.1) %>%
  filter(ER_2.5pct<0)%>%
  filter(K600_daily_Rhat<1.1)

mean_k_mod <- mean(met.clean$K600_daily_mean)
mean_k_mod

## [1] 21.11758

mean_k_measure <- mean(measured_K$K600)
mean_k_measure

## [1] 22.96813

## [1] -0.863

## [1] 0.185

Compared modeled and measured K600

Plots for (1) measured v modeled K600 and flow and (2) logK600 and log(flow+1).

Could be one poor measurement at the highest flow for measured gas exchange. But in general the modeled K600 does seem similar to the measured, which is kind of nice to see how robust the pool_K600 = "normal" is getting at K600.

mean_k_mod <- mean(met.clean$K600_daily_mean)
mean_k_mod

## [1] 21.11758

mean_k_measure <- mean(measured_K$K600)
mean_k_measure

## [1] 22.96813

k.q <- met.clean %>%
  left_join(mod.env.ag, by = "date") %>%
  full_join(measured_K, by = c("date", "source",
                                 "K600_daily_mean"= "K600",
                                "discharge"= "Q_cms")) %>%
  ggplot(aes(x = discharge, y = K600_daily_mean, col = source, shape=source))+
  geom_point( size = 2) +
  theme_bw()+
  #scale_x_log10(limits = c(20, 800))+
   scale_color_viridis_d(option = "viridis") 


k.q_log <- met.clean %>%
  left_join(mod.env.ag, by = "date") %>%
  full_join(measured_K, by = c("date", "source",
                                 "K600_daily_mean"= "K600",
                                "discharge"= "Q_cms")) %>%
  ggplot(aes(x = log(discharge+1), y = log(K600_daily_mean), col = source, shape=source))+
  geom_point( size = 2) +
  theme_bw()+
  #scale_x_log10(limits = c(20, 800))+
   scale_color_viridis_d(option = "viridis") 



plot_grid(k.q, k.q_log, ncol = 2)

possible_k <- (mean_k_mod+mean_k_measure)/2
possible_k

## [1] 22.04286

Full run with binned K600 informed by both measured and modeled priors

Here is the run configuration for full model:

mm_name(type = 'bayes', pool_K600 = "binned", err_obs_iid = TRUE, err_proc_iid = TRUE, ode_method = "trapezoid", deficit_src = 'DO_mod', engine = 'stan')

Quick plot of flow bins given flow range:

Where dashed vertical lines correspond to the prior locations for flow bins in bayes_specs_new$K600_lnQ_nodes_centers <- log_bins

Plots made on filtered data: met.clean filtered for days with GPP_daily_Rhat<1.05,ER_daily_Rhat<1.05, K600_daily_Rhat <1.05, as well as (GPP_97.5pct>0) and (ER_2.5pct<0).

met.clean <- met.full %>%
  filter(GPP_daily_Rhat<1.05)%>%
  filter(GPP_97.5pct>0)%>%
  filter(ER_daily_Rhat<1.05) %>%
  filter(ER_2.5pct<0)%>%
  filter(K600_daily_Rhat<1.05)

mean_k_mod <- mean(met.clean$K600_daily_mean)
mean_k_mod

## [1] 22.7534

mean_k_measure <- mean(measured_K$K600)
mean_k_measure

## [1] 22.96813

KER_cor <- round(cor(met.clean$ER_daily_mean, met.clean$K600_daily_mean, use = "complete.obs"),3)
print(KER_cor)

## [1] -0.817

KGPP_cor <-round(cor(met.clean$GPP_daily_mean, met.clean$K600_daily_mean, use = "complete.obs"),3)
print(KGPP_cor)

## [1] 0.347

The vertical dashed is the overall mean modeled K600 in the box plot.

GPP and ER

## [1] -0.734

Final thoughts.

The direction of the K600 ~ flow relationship looks more logical, where K600 increases with flow. ER and K600 are negatively correlated (-0.543), GPP and K600 are negatively correlated (-0.003). but, less strongly relative to the lower reach (GBL). The relationship between K600 and flow appears to be positive but still inflected in a strange way.

However I’m still think we should be cautious in over interpreting ER trends.

The mean modeled K600 and measured gas exchange are essentially the same 22.

Final cleaned model output:

Where GPP is in blue and ER is in orange, and the black points represent NEP.

Table summarizing excluded data from cleaned model output:

Metric	Number of Days	Percent of FULL DO days
Total DO days (FULL: >= 48 half-hour intervals)	860	100.0
Days with ANY metabolism output (on full DO days)	773	89.9
Days with NO metabolism output (on full DO days)	87	10.1
Days with CLEAN metabolism (GPP+ER+K600 all pass)	284	33.0
Days with CLEAN GPP (Rhat<=1.05 & GPP_2.5pct>=0)	382	44.4
Days with CLEAN ER (Rhat<=1.05 & ER_97.5pct<=0)	764	88.8
Days with CLEAN K600 (Rhat<=1.05)	642	74.7
Days with UNREASONABLE GPP (Rhat>1.05 & GPP_97.5pct<0)	0	0.0
Days with UNREASONABLE ER (Rhat>1.05 & ER_2.5pct>0)	0	0.0
Days with poor K600 (Rhat>1.05)	131	15.2

new_df_clean <- met_class %>%
  filter(clean_all==T)

KER_cor <- round(cor(new_df_clean$ER_daily_mean, new_df_clean$K600_daily_mean, use = "complete.obs"),3)
print(KER_cor)

## [1] -0.04

KGPP_cor <-round(cor(new_df_clean$GPP_daily_mean, new_df_clean$K600_daily_mean, use = "complete.obs"),3)
print(KGPP_cor)

## [1] -0.434

GPP_ER_cor <-round(cor(new_df_clean$GPP_daily_mean, (new_df_clean$ER_daily_mean*-1), use = "complete.obs"),3)
print(GPP_ER_cor)

## [1] 0.558

Session info

R version 4.4.2 (2024-10-31)

Platform: aarch64-apple-darwin20

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: plotly(v.4.10.4), kableExtra(v.1.4.0), knitr(v.1.50), streamMetabolizer(v.0.12.1), ggpubr(v.0.6.0), readxl(v.1.4.3), zoo(v.1.8-12), cowplot(v.1.1.3), viridis(v.0.6.5), viridisLite(v.0.4.2), dataRetrieval(v.2.7.17), lubridate(v.1.9.4), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.1.0), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.3.0), ggplot2(v.3.5.2) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): DBI(v.1.2.3), gridExtra(v.2.3), rlang(v.1.1.6), magrittr(v.2.0.3), e1071(v.1.7-16), compiler(v.4.4.2), mgcv(v.1.9-1), systemfonts(v.1.1.0), vctrs(v.0.6.5), pkgconfig(v.2.0.3), crayon(v.1.5.3), fastmap(v.1.2.0), backports(v.1.5.0), labeling(v.0.4.3), pander(v.0.6.5), deSolve(v.1.40), rmarkdown(v.2.29), tzdb(v.0.4.0), bit(v.4.5.0.1), xfun(v.0.53), cachem(v.1.1.0), jsonlite(v.2.0.0), broom(v.1.0.7), parallel(v.4.4.2), R6(v.2.6.1), bslib(v.0.9.0), stringi(v.1.8.7), RColorBrewer(v.1.1-3), car(v.3.1-3), jquerylib(v.0.1.4), cellranger(v.1.1.0), Rcpp(v.1.1.0), Matrix(v.1.7-1), splines(v.4.4.2), timechange(v.0.3.0), tidyselect(v.1.2.1), rstudioapi(v.0.17.1), dichromat(v.2.0-0.1), abind(v.1.4-8), yaml(v.2.3.10), lattice(v.0.22-6), plyr(v.1.8.9), withr(v.3.0.2), evaluate(v.1.0.5), rLakeAnalyzer(v.1.11.4.1), sf(v.1.0-21), units(v.0.8-7), proxy(v.0.4-27), xml2(v.1.4.0), pillar(v.1.11.0), carData(v.3.0-5), KernSmooth(v.2.23-24), generics(v.0.1.4), vroom(v.1.6.5), hms(v.1.1.3), scales(v.1.4.0), class(v.7.3-22), glue(v.1.8.0), lazyeval(v.0.2.2), tools(v.4.4.2), data.table(v.1.16.4), ggsignif(v.0.6.4), LakeMetabolizer(v.1.5.5), grid(v.4.4.2), crosstalk(v.1.2.2), nlme(v.3.1-166), Formula(v.1.2-5), cli(v.3.6.5), svglite(v.2.1.3), gtable(v.0.3.6), rstatix(v.0.7.2), sass(v.0.4.10), digest(v.0.6.37), classInt(v.0.4-11), htmlwidgets(v.1.6.4), farver(v.2.1.2), htmltools(v.0.5.8.1), lifecycle(v.1.0.4), httr(v.1.4.7), unitted(v.0.2.9) and bit64(v.4.5.2)